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<D . ABSTRACT 

[-^ I We tackle the problem of recovering an unknown signal observed in an ill-posed inverse 

problem framework. More precisely, we study a procedure commonly used in numerical 

E"' . analysis or image deblurring: minimizing an empirical loss function balanced by an penalty, 
OO 

! acting as a sparsity constraint. We prove that, by choosing a proper loss function, this 
> . 

^ ' estimation technique enables to build an adaptive estimator, in the sense that it converges 
at the optimal rate of convergence without prior knowledge of the regularity of the true 

, solution. 

§ : INTRODUCTION AND NOTATION 
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In this article we are interested in recovering an unobservable signal Xq based on obser- 
vations 

y{ti) = F{xo){U)+ei, (1) 



^ I where F : X Y is a. linear functional, with X, Y Hilbert spaces and ti,i = 1 



, n IS a 

fixed observation scheme. Xq : R ^ R is the unknown function to be recovered from the 
data y{ti), i = 1, . . . ,n. The regularity condition over the unknown parameter of interest is 
expressed through the assumption xq G X and will be made precise later in Section 3. We 
assume that the observations y{ti) G R and that the observation noise Ei are i.i.d. realizations 
of a certain random variable e. Throughout the paper, we shall denote y = 
We assume F is Frechet differentiable and ill-posed in the sense that our noise corrupted 
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observations might lead to large deviations when trying to estimate xq. 

In a deterministic framework, the statistical model ([T]) is formulated as the problem of 
approximating the solution of 

Fix) = y, 

when y is not known, and is only available through an approximation y^, 

\\y-y'\\^S. 

It is important to remark that whereas in this case consistency of the estimators depends on 
the approximation parameter 6, in ([1]) it depends on the number of observations n. 

The best approximation of Xq is = F'^y, where F'^ is the Moore-Penrose (general- 
ized) inverse of F. We will say the problem is ill-posed if is unbounded. This entails that 
F'^iy^) is not close to Hence, the inverse operator needs to be, in some sense, regularized. 
Regularization methods replace an ill-posed problem by a family of well-posed problems. 
Their solution, called regularized solutions, are used as approximations of the desired solu- 
tion of the inverse problem. These methods always involve some parameter measuring the 
closeness of the regularized and the original (unregularized) inverse problem. Rules (and 
algorithms) for the choice of these regularization parameters as well as convergence proper- 
ties of the regularized solutions are central points in the theory of these methods, since they 
allow to find the right balance between stability and accuracy. 

For this we consider penalized M-estimators minimizing quantities of the form 

Xn = argmin (7„(?/ - F{x){t)) + a„pen(a;, X)) , (2) 

where A" is a specific set, 7n(-) is an empirical loss-function, pen(., .) is a penalty over x in 
X, and On € 6 is a decreasing sequence all of which will be defined precisely later. The 
idea of penalized M-estimators is to find an estimator close enough to the data, close in 
the sense defined by 7„ and with a regularity property induced by the choice of the penalty 
pen. The smoothing sequence a„ balances the two terms. The greater a„, the smoother 
the estimator will be, while the smaller a„ the closer the estimator will be to the data, 
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maybe leading to a too rough estimate. Moreover the penalty should be chosen without 
prior knowledge of the regularity of the function to be estimated, Xq, in order to give rise to 
adaptive estimator. Adaptivity here has to be understood in the sense that the estimator 
converges at the optimal rate without knowing a priori its regularity. It is a theoretical 
alternative to cross validation methods developed for instance in Dey et al. (1996). We 
point out that in numerical analysis, adaptivity is not a key issue as in statistics. Indeed the 
smoothing sequence is often selected using such posterior techniques, see Tautenhahn and 
Jin (2003) or Kaltenbacher (2000). In a deterministic setting, the rates are not altered. This 
not the case when observing the data in a white noise framework, which justifies the need 
for fully adaptive methods. 

Different choices of penalty have been investigated in the literature. The traditional 
choice of a quadratic penalty defines the Tikhonov regularized estimator whose behaviour is 
well studied but which does not lead to adaptive estimation. For general references about 
this estimator, we refer to Tikhonov et al. (1998), Bissantz et al. (2004), or Engl (2000). A 
penalty on the number of non zero coefficients leads to hard-thresholded or projection esti- 
mators, whose asymptotic behaviour is studied in Kaltenbacher (2000), Mair and Ruymgaart 
(1996) or Engl et al. (1996). Adaptive estimators can yet be built with both methodologies, 
but by considering model selection techniques, as done in Loubes and Ludcna (2005). But 
adaptivity means, in that case, that the estimator behaves as well as the best estimate ob- 
tained over a fixed class of estimators, i.e a collection of models. This property is generally 
expressed through an oracle inequality. Hence, there might be a bias if the true solution is 
not well approximated by the sieves. 

In this article, we tackle the problem of the asymptotic behaviour of the estimator ob- 
tained with a penalty. Indeed, over the last decade, penalty has been more and more 
used in a large variety of fields. Indeed such a penalty selects sparse signals in a smoother 
way than hard-thresholding penalty and can also be easily implemented. Contrary to differ- 
entiablc penalties for which adaptivity implies selecting the smoothing sequence among a set 
of possible choices, there is an optimal choice of the trade-off parameter when using a soft- 
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thresholding penalty. And this optimal choice does not depend on the unknown regularity 
of the parameter of interest, which enables adaptive estimation. Hence norm penalty is 
used in estimation with soft-thresholding estimators in Loubes and van de Geer (2002) or 
Loubes (2007), and in inverse problems in Daubechies et al. (2004) or Cohen et al. (2003). 
We point out that in Cohen et al. (2003), the estimator is the soft thresholded version of the 
estimation which is used, and the properties of the estimator are studied using a sequential 
version of model ([1]) together with a projection method. 

In this article, we construct a penalized estimator with a penalty, with an appropriate 
loss function depending on the operator. We show that such an estimator converges and is 
adaptive over a class of Besov spaces for a certain class of ill-posed problems. When trying 
to minimize a standard empirical quadratic loss function together with a softhresholding 
penalty, as often done in numerical analysis literature, the results are different. The estimator 
may be inconsistent and its rate of convergence is rather slow. We also provide the rate of 
convergence in this case and stress the advantages of choosing a loss function depending on 
the operator to build a more efficient estimator. 

The article falls into the following parts. Section 1 presents the model and the overall 
assumptions. The estimation procedure and its general efficiency are described in Section 
2. Rates of convergence and adaptivity under smoothness assumptions are also discussed in 
Section 3. The section 4 is devoted to the analysis of the performances of the usual penalized 
least squares estimator for inverse problems. Simulations are conducted in Section 5. 

1. INVERSE PROBLEM MODEL 

Consider the following inverse model: 



where F : X ^ Y is a known linear operator, whose adjoint will be denoted F*. Set S the 
Dirac function and define Qn the empirical measure of the covariables as Qn = ^ J2i=i ■ 
Throughout all the paper, the estimation errors will be given with respect to the L2{Qn)- 




(3) 
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norm defined, for all functions y G Y, by 

1 " 

n ^-^ 

1=1 

The corresponding empirical scalar product is given by < y,e >n— \ Yll=i ^iVi^i)- 

As often F is not of full rank, so the singular value decomposition (SVD) of the operator 
S is then a useful tool. 

Let iXj','fij,'4^j)j^i be a singular system for the linear operator F, that is, Fipj — Xj(pj and 
F*(pj = Xji/jj] where {A^j^^i are the non zero eigenvalues of the selfadjoint operator F*F, 
considered in decreasing order. Furthermore, {■(/'j}j=i,...,n and {(Pj}j=i,...,n ai'e a correspond- 
ing complete orthonormal system with respect to of eigenvectors of F*F and FF*, 
respectively. For general linear operators with an SVD decomposition, we can write for all 
(x,y)eXxY 

n 

Fx^^Xj{x,ipj)^j (4) 

n 

For y in the domain of F\ 'D{F^), the best-approximate solution has the expression 

Note that for large j, the term 1/Xj grows to infinity. Thus, the high frequency errors 
are strongly amplified. This amplification measures the difficulty of the inverse problem, the 
faster the decay of the eigenvalues, the more difficult is the inverse problem. In this paper 
we will tackle the problem of polynomial decay of eigenvalues, which can be described by 
the following assumption 

Index of ill-posedness Assume that there exists an index i > 0, called the index of ill- 
posedness of the operator F, following notations in Engl et al. (1996), such that 
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This difference with standard regression model for which t = yields other optimal rates 
which are usual in statistics. In Section 3, we will compare the rates obtained by our 
estimator to these optimal rates of convergence. 

A penalized M-estimator is defined using a distance d, between the observations y = 
{y(ti), . . . ,y{tn)) G Y and a function x E X, and a penalty. Hence we shall study an 
estimator of the following type 



Xn = argmm 



1 " 

-J2d{y{U),x{U)) + I^( 



n 
1=1 



(6) 



where /i is a smoothing sequence and /^(.) is taken to be the soft-thresholding penalty. 
Indeed for = (/ij), j = 1, . . . , and Vx = X]j=i ^j'^i ^ -^i -^Ai(^) ■= TTj=i l/^j^il 
weighted norm of the function x. In the direct case where F = Id, and for a quadratic loss 
function, i.e 



^ n 1 

- V d{y{U),x{U)) = - V \y{U) - x{U) I 
n ^ — ^ n ^-^ 

i=l i=l 
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the estimator can be computed explicitly and is called the soft-thresholded estimator, as 
pointed out in Loubes and van de Geer (2002). 

In the inverse problem literature, such an estimator is often used in image recognition or 
in geophysics. Indeed the norm is well fitted to handle such signals. In such cases, it is 
not possible to solve the corresponding minimization issue and the asymptotic properties of 
the estimator are not known. Hence, in the following section, we propose an adequate choice 
of loss function which enables to build a sparse estimator converging at an optimal rate of 
convergence. Moreover this estimator is adaptive, thanks to the sparsity property of the 
norm. We point out that adaptivity means that the regularity of the function xq is unknown 
but that the estimator still achieves the optimal rate of convergence. Nevertheless, the 
operator F is assumed to be known, as well as the degree of ill-posedness. This assumption 
is common in the statistical literature on inverse problems. 

2. ESTIMATION USING SOFT THRESHOLDING PENALTY 

In this section, we investigate the classical inverse regression model ([1]) with independent 
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errors ei,...,e„ with zero expectation and finite variance . Assume moreover that the 
following condition over the observations errors holds 

Error bound (A): Suppose that for some constant K < oo, the errors ei, . . . , e„ satisfy 

max Eexp[e,V^1 ^ ^- 

i=l,...,n 

For a choice of penalty /^(x) = Yl]=i hI^jI^ be the penalized estimator defined as 



arg mm 



E 

.i=i 



2 

<y-F{x),^>^ 



(7) 



with = (/Xj, J = 1, . . . , n) a sequence of real numbers. 

This penalized estimator mimics the soft-thresholdcd estimator, as pointed out in Loubes 
and van de Geer (2002). Contrary to the direct case where the thresholding level can 
be chosen equal to a constant = in this case we consider a threshold that changes 
at each reconstruction level j, /i — {/ij) and which depends on the nature of the inverse 
problem. Indeed, selecting the coefficients without considering the effect of the inverse 



problem is too rough and the usual choice in nonparamctric estimation fi = 2cy^^ may 
lead to inconsistent estimator. Hence, in this this work, define the smoothing sequence as 
Vj = 1, . . . , n jij :— 2^-^^^", for c a given constant. 

Write Xq — YTj=i ^jfii^j consider = Yl^=i ^j,*'^3-i hard-threshold version of the 
true function xq, defined as 



1 ^j,o> 




> 






if ki,o| 







"^3:* — \ 1 j — 1 , . . . , 77/. 



We will now establish an upper bound for — Xo||n which depends on the performance of 

the oracle 11 II This will enable us to get rates of convergence for ill-posed inverse 

problems and to prove adapt ivity of the estimation procedure. 
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In the theorem we write 



1 " 

Theorem 1. Let Bn he the set 



i=l 



I log 71 

Bn = { max \Vj\ ^ c\ }. 

j=l,...,n \ n 

Consider the set of indexes Jn = Cardjj, \xjfi\ > fij}. Then on we have 

'logrz 1% ~ ^j,o| 



Proof. First note that the empirical contrast can be rewritten in a different way for all x G X 



E 



A, 



<y- F{x),ipj >, 



1 



— < e, Lpj >n 
^3 



2 n 

=5: 

2 n 



[< y - F(Xo), >n + < F{xo) - F{x), >„] 



< e, V^j >n 



Set V^- =< e, yjj >„= ^ X]r=i ^»V^i(^i)- Using Definition ([7j) and previous remark, it implies 
that 



E 



< ?/ - F{Xn),(Pj >r 



<y - F{x^),ipj >n 



I X j X j ^ 



+ I{nx^) 



^ \\x 



,-Xo||^ + 2( maxjVSlj^ ^ 



Now use the following decomposition 



J h 



to get that on the set 



W^n - XoWl + I^,{Xn) ^ \\X^ - Xq\\1 + 2cy ^ ^ + I^,{x^) 

j=l ■? 

II- 112 I- I ^ II 112 ^ jhgn -s-^ \xj — Xj J ^ /logn IxJ + b,- * 

+ ~ l%l) + 5^ A*jl^j,*l- 

Point out that 

{Xjfi, if j e Jn 
0, ifj^Jn 

Then, for a choice /jLj — 2i2-^log n/n and using \ \xj\ — \xj^o\ \ ^ |xj — Xj^o\, we obtain that: 



logn 


E 




n 




A, 






logn 






n 


E 


A, 







ll^n - a^oll^ ^ Ik* - XoWl + 4c 

which proves the result. □ 

Corollary 1. Under assumption (A), it follows from e.g. van de Geer (2000) , Lemma 

8.2), that for a constant c depending on K, 



o I IT/I ^ \ ^ r logn . . 

P max \VA > c\ ^ cexp —\. (9) 

' j=l,...,n ' \ n I & 

Thus, we obtain for two positive finite constants Ci and C2 

logn. 



( 2 ii2 logn ^-^ 1 \ 

P I \\Xn - Xq\\^ > Ci\\x^ - XoWn + ^2^^ J2 ) ^ Cexp[ 



c2 



We point out that previous bound is as sharp as the equivalent one for direct estimation 
problems, see for instance Loubes and van de Geer (2002). The choice of the smoothing 
sequence /i does not depend on the regularity of the unknown function xq. Hence we expect 
adaptivity under regularity conditions. It still depends on the distribution of the errors, 
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since the constant c in ([9]) may be large. As a consequence, if the errors have heavy tails, 
the rate of convergence of the penalized estimator may be slow. 

Nevertheless, this estimator works with real observation and can handle a large variety 
of observation noise since Assumptions (A) are rather weak. 

Proof. First note that Cauchy-Schwartz inequality entails that, on the set Bn, we have 



I- \\2 ^ A /log":, 



n 



Since for all 1 > 7 > 0, we get 2xy ^ -x + jy , we obtain that 

(i-7)ll*„-x„ll;<lk.-:^olli; + ^^i^E4 

leading to ([1]) for proper constant choices. 



□ 



3. RATES OF CONVERGENCE UNDER SOURCE CONDITION 

In this section, we illustrate the consequences of Corollary [1] for functions xq belonging 
to some special smoothness sets. Consider the set of functions defined by two parameters, a 
smoothness parameter s and a moment parameter < p < 2 as 

n n 

Xs,p = {x = ^ 1}. (10) 

i=i i=i 

Such sets are balls of Besov bodies associated to the Besov spaces i?pp([0, 1]). These spaces 
are intrinsically connected to the analysis of curves since the scale of Besov spaces yields the 
opportunity to describe the regularity of functions, with more accuracy than the classical 
Holder scale. General references about Besov spaces are Besov et al. (1978). Consider a 
wavelet basis {'ipjk)j,k with regularity r such that r ^ s. Then a Besov norm is equivalent to 
an appropriate norm in the sequence space, that is the space of the wavelet coefficients. If 
Xjk are the wavelet coefficients of a function x = 'Ylj hence the ball with radius 1 

of Besov space i?pp([0, 1]) can be fully characterized by the Besov semi-norm 



^+00 






V 




E 




("£■■■; 


1 












I 



^ 1. 
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as proved in Hardle et al. (1998). In the Besov space interpretation Xg^p with n = oo 
corresponds (in the sense of norm equivalence) to a Besov ball in the space Bpp{[0, 1]). 

Here consider the special case where p is such that ^ = | + 2t+T- This choice corresponds 
to the set for which classical soft thresholded estimators for inverse problems are optimal, 
see for instance Theorem 3.2 in Cohen et al. (2003). For that choice of parameters, we have 
jp{s+i/2~i/p) _ j^-2sp/(2t+i)^ Hence previous set can then be rewritten as 

n n p 



Xs,p = {X = Y^Xji)j, ^ 1}. 



. , • 1 \ 2t+l 

Such regularity condition can be interpretated as a source set condition, used in the literature 
of deterministic inverse problems, see for instance Engl et al. (1996), Darolles et al. (2003) or 
Fermin et al. (2005). Such spaces link the decay of the Xj's, the coefficients of the unknown 
function in the SVD basis with the decay of the X/s, the eigenvalues of the operator. 

The following theorem gives the rate of convergence of the penalized estimator for 
inverse problems. 

Theorem 2. Assume that there are s and < p < 2 such that xq G Xg^p, with ^ = ^ + 27+T- 
Then we get the following estimation error 



|£n-Xo||^ = Op| (^J j. (12) 

Proof. Starting from Corollary [T], we have to bound the two terms || II „, which stands 

for a bias term, and y^,v--7 tt which stands for a variance term. For this we will take 
Jn = {j) Xjf) ^ make use of the assumption X]j>i A~^'*^^^^*'*""'^'*x^q ^ 1. 
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First write 



I ||2 \ ^ 2 



2tps 2tp. 



^x'jj2t+ij 2t+i 



[ogn\ 2 „ .2lp± .t(2-v)~2lPi. 



n / \ n 



using repeatdly ^ = 1 + iFfT ^^'^ ^(^ - P) - ^ = 0- 
On the other hand we get 



\ 2 



■t(2-p) -tp 

jeJn ^ jeJn 



E 

< ( — 1 y^x^j^j 



logn 



*{2-p)-iff 



n 



logn^ 

Previous upper bounds lead to (fT2|) . concluding the proof. □ 



We point out that we obtain the optimal rate of convergence for ill-posed inverse problems, 
i.e (n) 2s+2t+i ^ up to a logarithmic factor. Hence, the estimator (jTj) is adaptive with respect 
to the parameter s within the range of Besov spaces Bpp{[0, 1]) with ^ = | + 2t+T 
empirical quadratic loss. This result is the same as the one obtained in Cavalier et al. (2002), 
when working in the sequential model. Hence we provide a new estimator which achieves 
optimal rates of convergence and which can be easily computed, as shown in Section 5. 

4. COMMENTS ON PENALIZED LSE FOR INVERSE PROBLEMS 

The estimator obtained in ([7j) uses a specific loss function adapted to the particular ill- 

2 

posed problem, namely dn{y,x) := J2j=i < V ~ ^{x), ^ >„ . This loss function leads to 
optimal rate of converge but its main drawback is that the knowledge of the operator F 
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and its SVD {Xj;(fj,ipj), j ^ 1 are needed. In numerical analysis for inverse problems, a 
classical deblurring procedure involves minimizing the usual quadratic empirical loss function 
together with the penalty, see for instance Daubechies et al. (2004) or Cohen et al. (2003). 
More precisely an estimator is defined as 

x„ = arg min {\\y - F{x)\\l + If,{x)) . (13) 

This estimation procedure is used and its weak consistency is well-known, even if it is known 
to provided sometimes an inconsistent estimate. However this estimator can be easily imple- 
mented in an iterative procedure, similar to a gradient descent algorithm with a data driven 
step. Such algorithm is widely used in image deblurring for instance. We refer to Daubechies 
et al. (2004) for more references. 

The following theorem gives conditions to ensure consistency and provides rates of con- 
vergence for the empirical loss function. 

Theorem 3. Assume that the errors e satisfy to the condition (A). Assume that there exists 
a roughness parameter < p < 2 such that 

n 

J2\^,,o\'<l. (14) 

Hence the estimator defined in f|T3|) is consistent as soon as 1 — p/2 — 2t > and converges 
at the following rate of convergence 



Xn-Xo\\l = Op{\og'-Hn){-] 1. (15) 



n 



l-f-2t^ 



Proof. Using the definition of the estimator, first point out that on the set -B„ = {maxj=i \ Vj\ 



cy^^^}, we have following the guidelines of the proof of Theorem [T] 

\\y - F{x^)\\l + /^(5„) ^ \\y - F{x^)\\l + /^(xj, 
leading for the oracle and a choice 



Jn = {j, \xj,o\ ^ c^y-^^}, 
Xj \ n 

13 



2 



a Jr. 

where \Jn\ stands for the cardinal of the set Jn- Now using that 

J 3&Jn J&Jn 



we get the following bound \/\Jr\ ^ (k^l^i)'''^'^- Finally, using that 2xy ^ + for all 
positive 7, we obtain that on the set i?„ we get 

i-§-2r 



log^ 2 (72 



1\- 2 



Assumption (A) ensures that P(-B^) ^ 0, which concludes the proof. □ 

We point out that the estimator x„ is convergent when the index of ill-posedness is such 
that t < \ — \- This implies that for t > 1/2, i.e for severe ill-posed problems, the estimation 
procedure does not lead to consistent estimates. Even in the case where the estimator 
converges, its rate of convergence is less than the rate of the optimal estimator x„. This 
difference comes from the fact that the standard quadratic loss function is not well fitted to 
handle inverse problems since the extra term in A| = entails a loss of order That 
is the reason why we propose a loss function in the definition ([7]) which one the one hand 
gets rid of this term and other hand induces a bias error which can be balanced by the 
penalty, which enables to achieve optimality. 

We can argue that the regularity assumption fll4p is more general than Assumption fllOp . 
However, increasing the regularity of the function to be estimated does not help increasing 
the corresponding rate of convergence of the estimator x^. Indeed the limiting factor in 
n~^* comes from the rate of decay of the eigenvalues and not the decay of the coefficients 
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Figure 1: Index of Ill-posedness t = .2 

Sj^Jn ^o,j- Hence estimator a:„ is outperformed by estimator x„. 

However Theorem [3] gives the range of apphcation of the standard penahzed estimator 
and provides a better understanding of its behaviour. 

5. NUMERICAL RESULTS 

In this section, we apply our estimation procedure to simulated data obtained using a 
sequence model. The function we wish to reconstruct is t —>■ sin(l/(.l + 1), observed in an 
ill-posed settings with t standing for the index of ill-posedness. We take n = 200 observations 
and the observation noise is a Gaussian white noise with snr = 2. 

We plot in straight lines the true signal to be recovered, in dotted lines the estimator while 
the observations are represented by crosses. We will consider two cases depending on the 
ill-posedness of the operator, with an easy inverse problem t < | in Figure [1] and a severe 
ill-posed problem t > ^ in Figure [H The two figures are significant realizations of 50 
replications. 
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50 100 150 200 250 

Figure 2: Index of Ill-posedness t = 1.5 

We can see that in both cases the estimator x„ is consistent and provides a quite good ap- 
proximation of the unknown function. Hence the estimation procedure given in (JTj) provides 
a good estimator of ill-posed inverse problems when nothing is known about the regularity 
of the function. Its main drawback is that a very good knowledge of the operator is needed, 
but this is the case in most of the denoising procedures for inverse problems. 
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